t 2 g $$ {t}_{2g} $$ d $$ d $$ orbital ordering patterns in KBF 3 $$ {}_3 $$ (B = Sc, Ti, Fe, Co) perovskites.

The orbital ordering (OO) resulting from the partial occupancy of the t 2 g $$ {t}_{2g} $$ d $$ d $$ subshell of the transition metals in KBF 3 $$ {}_3 $$ (B = Sc, Ti, Ffe, Co) perovskites, and the many possible patterns arising from the coupling between the B sites, have been investigated at the quantum mechanical level ( all electron $$ \mathrm{all}\ \mathrm{electron} $$ Gaussian type basis set, B3LYP hybrid functional) in a 40 atoms supercell. The numerous patterns are distributed into 162 classes of equivalent configurations. For each fluoroperovskite, one representative per class has been calculated. The four compounds behave similarly: an identical dependence of the energy and volume (or cell parameters) on the OO pattern; the spanned energy interval is small (1 to 2 mE h $$ {}_h $$ per formula unit), suggesting that most of the configurations are occupied at room and even at low temperature. A linear model, taking into account the relative orbital order in contiguous sites, reproduces the energy order in the full set for each compound, suggesting that it could be used for studying OO in larger supercells.

The d Orbital Multi Pattern Occupancy in a Partially Filled d Shell: The KFeF3 Perovskite as a Test Case.

The occupancy of the d shell in KFeF3 is t2g4eg2, with five α and one ß electrons. The Jahn-Teller lift of degeneracy in the t2g sub-shell produces a tetragonal relaxation of the unit cell (4.09 vs. 4.22 Å, B3LYP result) not observed experimentally. In order to understand the origin of this apparent contradiction, we explored, with a 2 × 2 × 2 supercell (40 atoms per cell), all possible local structures in which contiguous Fe atoms have a different occupancy of the t2g orbitals with the minority spin electron. A total of 6561 configurations (with occupancies from (8,0,0) to (3,2,2) of the 3 t2g orbitals of the 8 Fe atoms) have been explored, with energies in many cases lower (by up to 1550 µEh per 2 Fe atoms) than the one of the fully ordered case, both for the ferromagnetic and the anti-ferromagnetic solutions. The results confirm that the orientation of the ß d electron of Fe influences the electrostatics (more efficient relative orientation of the Fe quadrupoles of the d shell) of the system, but not the magnetic interactions. Three hybrid functionals, B3LYP, PBE0, and HSE06, provide very similar results.

Strategies for the optimization of the structure of crystalline compounds.

When different proposals exist (or can reasonably be formulated) for the size of the unit cell (in terms of number of atoms) and space group of crystalline compounds, a strategy for exploring with simulation methods the various cases and for investigating their relative stability must be defined. The optimization schemes of periodic quantum mechanical codes work in fact at fixed space group and number of atoms per unit cell, so that only the fractional coordinates of the atoms and the lattice parameters are optimized. A strategy is here presented, based on four standard tools, used synergistically and in sequence: (1) the optimization of inner coordinates and unit cell parameters; (2) the calculation of the vibrational frequencies not only at Γ , but also at a set of k → points (in the example presented here they are eight, generated by a shrinking factor 2), looking for possible negative wavenumbers. The latter correspond to maxima, rather than minima, along the coordinate described by the corresponding normal mode; (3) the exploration of the total energy along the mode with negative wavenumber, looking for the minimum of the curve; (4) the identification of the new space group corresponding to the reduced symmetry resulting from the previous step. The strategy is illustrated with reference to the KMnF3 perovskite compound, for which many space groups are proposed in the literature, ranging from cubic Pm 3 ¯ m to tetragonal P 4 m bm or I 4 m cm and orthorhombic (Pnma and Cmcm) down to monoclinic (P21 /m). The corresponding primitive cells contain 5, 10, and 20 atoms in the various cases, and the point symmetry reduces from 48 to 4 operators. In nature, the KMnF3 phase transitions also include the magnetic phases. For simplicity, here we limit the analysis to the ones that take place between ferromagnetic phases, as they are sufficiently rich for illustrating the proposed strategy. As the total energy differences involved can be as small as, say, 10-50 µHartree, a high numerical accuracy at each one of the steps mentioned above is required. The present calculations, performed with the CRYSTAL code, by using an all electron basis set and the Hartree-Fock and B3LYP functionals, document such an accuracy. The energy difference between the tetragonal I 4 m cm and cubic Pm 3 ¯ m phases is 225 µHartree, with a volume reduction of 0.58 Å3 ; the differences between the orthorhombic and tetragonal phases are an order of magnitude smaller, being 23 µHartree and 0.06 Å3 for total energy and cell volume, respectively.

The ferromagnetic and anti-ferromagnetic phases (cubic, tetragonal, orthorhombic) of KMnF3. A quantum mechanical investigation.

Many space groups are proposed in the literature for the KMnF3 perovskite (see, for example, Knight et al., J. Alloys Compd., 2020, 842, 155935), ranging from cubic (C) (Pm3̄m) to tetragonal (T) ( or I4/m) down to orthorhombic (O) (Pbnm). The relative stability ΔE of these phases, both ferromagnetic (FM) and antiferromagnetic (AFM), has been investigated quantum mechanically by using both the B3LYP hybrid functional and the Hartree-Fock Hamiltonian, an all-electron Gaussian type basis set and the CRYSTAL code. The O phase is slightly more stable than the T phase which in turn is more stable than the C phase, in agreement with experimental evidence. The C to T to O transition is accompanied by a volume reduction. The mechanism of stabilization of the AFM solution with respect to the FM one is discussed. Spin density maps and profiles, Mulliken charges, magnetic moments and bond population data are used for supporting the proposed mechanism. The IR and Raman spectra of the FM and AFM C, T and O cells are discussed; the only noticeable difference between the C, T and O spectra appears at wavenumbers lower than 150 cm-1. The effect of pressure is also explored in the 0-20 GPa interval. The stability order (O > T > C) at 0 GPa persists also at high pressure, and the differences between the phases increase.

The superexchange mechanism in crystalline compounds. The case of KMF3(M = Mn, Fe, Co, Ni) perovskites.

The ferromagnetic and antiferromagnetic wavefunctions of four KMF3(M = Mn, Fe, Co and Ni) perovskites have been obtained quantum-mechanically with the CRYSTAL code, by using the Hartree-Fock (HF) Hamiltonian and three flavours of DFT (PBE, B3LYP and PBE0) and anall-electronGaussian type basis set. In the Fe and Co cases, with d6and d7occupation, the Jahn-Teller distortion of the cubic cell is as large as 0.12 Å. Various features of the superexchange interaction energies (SIE), namely additivity, dependence on the M-M distance, on theMFM̂angle, and on the adopted functional, are explored. The contribution to SIE by the Coulomb, exchange and kinetic energy terms is analyzed. It is shown that, when using density functionals, SIE clearly correlates with the amount of exact (Hartree-Fock) exchange in the functional. The effect of SIE on the equilibrium geometry and volume of the unit cell is discussed, and it is shown that the key quantity is the spin polarization of the (closed shell) F ions along the M-F-M path. The effect of thismagneticpressureis evaluated quantitatively for the first time. The superexchange coupling constantJ, evaluated at the HF level and through the Ising model, underestimates the experimental values by about 60%-70%. The more sophisticated Yamaguchi model (that takes into account the contamination of the FM and AFM spin states) does not reduce the discrepancy. The B3LYP hybrid functional overestimates the experiments. These last are bracketed by HF and PBE0. For PBE, the overestimation is huge. Finally, Mulliken population data, charge and spin density maps and density of states are used to illustrate the electronic structure.

Interstitial carbon defects in silicon. A quantum mechanical characterization through the infrared and Raman spectra.

The Infrared (IR) and Raman spectra of various interstitial carbon defects in silicon are computed at the quantum mechanical level by using an all electron Gaussian type basis set, the hybrid B3LYP functional and the supercell approach, as implemented in the CRYSTAL code (Dovesi et al. J. Chem. Phys. 2020, 152, 204111). The list includes two 〈100〉 split interstitial IXY defects, namely ICC and ICSi , a couple of related defects that we indicate as IX IY , the so called C i C s 0 in its A and B form, as well as SiCi Si and Cs Ci Cs , in which the interstitial carbon atom is twofold coordinated. The second undergoes a large relaxation, and the final configuration is close to ICC Cs . Geometries, relative stabilities, electronic, and vibrational properties are analysed. All these defects show characteristic features in their IR spectrum (above 730 cm- 1 ), whereas the Raman spectrum is dominated, in most of the cases, by the pristine silicon peak at 530 cm-1 , that hides the defect peaks.

The CRYSTAL code, 1976-2020 and beyond, a long story.

CRYSTAL is a periodic ab initio code that uses a Gaussian-type basis set to express crystalline orbitals (i.e., Bloch functions). The use of atom-centered basis functions allows treating 3D (crystals), 2D (slabs), 1D (polymers), and 0D (molecules) systems on the same grounds. In turn, all-electron calculations are inherently permitted along with pseudopotential strategies. A variety of density functionals are implemented, including global and range-separated hybrids of various natures and, as an extreme case, Hartree-Fock (HF). The cost for HF or hybrids is only about 3-5 times higher than when using the local density approximation or the generalized gradient approximation. Symmetry is fully exploited at all steps of the calculation. Many tools are available to modify the structure as given in input and simplify the construction of complicated objects, such as slabs, nanotubes, molecules, and clusters. Many tensorial properties can be evaluated by using a single input keyword: elastic, piezoelectric, photoelastic, dielectric, first and second hyperpolarizabilities, etc. The calculation of infrared and Raman spectra is available, and the intensities are computed analytically. Automated tools are available for the generation of the relevant configurations of solid solutions and/or disordered systems. Three versions of the code exist: serial, parallel, and massive-parallel. In the second one, the most relevant matrices are duplicated on each core, whereas in the third one, the Fock matrix is distributed for diagonalization. All the relevant vectors are dynamically allocated and deallocated after use, making the code very agile. CRYSTAL can be used efficiently on high performance computing machines up to thousands of cores.

Substitutional carbon defects in silicon: A quantum mechanical characterization through the infrared and Raman spectra.

The infrared (IR) and Raman spectra of eight substitutional carbon defects in silicon are computed at the quantum mechanical level by using a periodic supercell approach based on hybrid functionals, an all electron Gaussian type basis set and the CRYSTAL code. The single substitutional C s case and its combination with a vacancy (C s V and C s SiV) are considered first. The progressive saturation of the four bonds of a Si atom with C is then examined. The last set of defects consists of a chain of adjacent carbon atoms C s i , with i = 1-3. The simple substitutional case, C s , is the common first member of the three sets. All these defects show important, very characteristic features in their IR spectrum. One or two C related peaks dominate the spectra: at 596 cm-1 for C s (and C s SiV, the second neighbor vacancy is not shifting the C s peak), at 705 and 716 cm-1 for C s V, at 537 cm-1 for C s 2 and C s 3 (with additional peaks at 522, 655 and 689 for the latter only), at 607 and 624 cm-1 , 601 and 643 cm-1 , and 629 cm-1 for SiC s 2 , SiC s 3 , and SiC s 4 , respectively. Comparison with experiment allows to attribute many observed peaks to one of the C substitutional defects. Observed peaks above 720 cm-1 must be attributed to interstitial C or more complicated defects.

On the Models for the Investigation of Charged Defects in Solids: The Case of the VN- Defect in Diamond.

Local charged defects in periodic systems are usually investigated by adopting the supercell charge compensated (CC) model, which consists of two main ingredients: (i) the periodic supercell, hopefully large enough to reduce to negligible values the interaction among defects belonging to different cells; (ii) a background of uniform compensating charge that restores the neutrality of the supercell and then avoids the "Coulomb catastrophe". Here, an alternative approach is proposed and compared to CC, the double defect (DD) model, in which another point defect is introduced in the supercell that provides (or accept) the electron to be transferred (subtracted) to the defect of interest. The DD model requires obviously a (much) larger supercell than CC, and the effect of the relative position of the two defects must be explored. A third possible option, the cluster approach, is not discussed here. The two models have been compared with reference to the VN- defect; for DD, the positive compensating charge is provided by a P atom. Three cubic supercells of increasing size (containing 216, 512, and 1000 atoms) and up to eight relative VN--P+ defect-defect positions have been considered. The comparison extends to the equilibrium geometry around the defect, band structure, charge and spin distribution, IR and Raman vibrational spectra, and electron paramagnetic resonance constants. It turns out that the CC and DD models provide very similar results for all of these properties, in particular when the P+ compensating defect is not too close to VN-.

Vibrational spectroscopy of hydrogens in diamond: a quantum mechanical treatment.

The electronic and vibrational features of the VHn (n = 1 to 4) family of defects in diamond (hydrogen atoms saturating the dangling bonds of the atoms surrounding a vacancy) are investigated at the quantum mechanical level by using the periodic supercell approach, an all electron Gaussian type basis set, hybrid functionals, and the Crystal code. Most of the results have been collected for supercells containing 64 atoms; however, in order to explore the effect of the defect concentration on both the IR and Raman spectra, supercells containing 216, 512 and 1000 atoms have also been considered in the VH4 case. For each system, all the possible spin states are considered; their relative stability, band structure, charge and spin density distributions are thoroughly described. All the investigated systems present specific IR and Raman spectra, with vibrational spectroscopic features that can in principle be used as fingerprints for their characterization. This is particularly true for the C-H stretching, that ranges between 2500 and 4400 cm-1. The stretching modes are strongly affected by anharmonicity that has been evaluated in this work; it turns out to be extremely sensitive to the H load and spin state of the system, and ranges from -335 cm-1 for VH1 to +85 cm-1 for VH4. All of the investigated defects have very low C-H stretching IR intensity, so that they essentially appear as silent, the exception being VH1. The situation is different for the Raman spectra: the stretching modes of all defects do have similar large intensity; unfortunately here it is the experimental evidence that is lacking.

Katoite under pressure: an ab initio investigation of its structural, elastic and vibrational properties sheds light on the phase transition.

The evolution under pressures up to 65 GPa of structural, elastic and vibrational properties of the katoite hydrogarnet, Ca3Al2(OH)12, is investigated with an ab initio simulation performed at the B3LYP level of theory, by using all-electron basis sets with the Crystal periodic program. The high-symmetry Ia3d phase of katoite, stable under ambient conditions, is shown to be destabilized, as pressure increases, by interactions involving hydrogen atoms and their neighbors which weaken the hydrogen bonding network of the structure. The corresponding thermodynamical instability is revealed by anomalous deviations from regularity of its elastic constants and by numerous imaginary phonon frequencies, up to 50 GPa. Interestingly, as pressure is further increased above 50 GPa, the Ia3d structure is shown to become stable again (all positive phonon frequencies and regular elastic constants). However, present calculations suggest that, above about 15 GPa and up to at least 65 GPa, a phase of I4[combining macron]3d symmetry (a non-centrosymmetric subgroup of Ia3d) becomes more stable than the Ia3d one, being characterized by strengthened hydrogen bonds. At low-pressures (between about 5 GPa and 15 GPa), both phases show some instabilities (more so for I4[combining macron]3d than for Ia3d), thus suggesting either the existence of a third phase or a possible phase transition of second order.

Elasticity of grossular-andradite solid solution: an ab initio investigation.

Grossular and andradite are garnet end-members stable under upper mantle conditions. We perform ab initio simulations to investigate the dependence of the bulk modulus on chemical composition of the grossular-andradite solid solution, Ca3Fe(2-2x)Al(2x)(SiO4)3. All-electron local basis sets of Gaussian-type orbitals and the hybrid B3LYP density functional are used. Our calculations predict a linear modulus-composition trend, in contrast to previous conjectures based on "heterogeneous" experimental measurements. We estimate the largest deviation from linearity to be about 0.5 GPa under ambient conditions, and to progressively reduce to less than 0.2 GPa at pressure P = 20 GPa. The bulk modulus is computed over the whole composition range 0 ≤x≤ 1 following two independent approaches: fitting energy-volume data to an equation-of-state and calculating elastic tensors. Results from the two methods are in perfect agreement, assuring consistency and high numerical accuracy of the adopted algorithms.

Symmetry and random sampling of symmetry independent configurations for the simulation of disordered solids.

A symmetry-adapted algorithm producing uniformly at random the set of symmetry independent configurations (SICs) in disordered crystalline systems or solid solutions is presented here. Starting from Pólya's formula, the role of the conjugacy classes of the symmetry group in uniform random sampling is shown. SICs can be obtained for all the possible compositions or for a chosen one, and symmetry constraints can be applied. The approach yields the multiplicity of the SICs and allows us to operate configurational statistics in the reduced space of the SICs. The present low-memory demanding implementation is briefly sketched. The probability of finding a given SIC or a subset of SICs is discussed as a function of the number of draws and their precise estimate is given. The method is illustrated by application to a binary series of carbonates and to the binary spinel solid solution Mg(Al,Fe)2O4.

On the use of symmetry in configurational analysis for the simulation of disordered solids.

The starting point for a quantum mechanical investigation of disordered systems usually implies calculations on a limited subset of configurations, generated by defining either the composition of interest or a set of compositions ranging from one end member to another, within an appropriate supercell of the primitive cell of the pure compound. The way in which symmetry can be used in the identification of symmetry independent configurations (SICs) is discussed here. First, Pólya's enumeration theory is adopted to determine the number of SICs, in the case of both varying and fixed composition, for colors numbering two or higher. Then, De Bruijn's generalization is presented, which allows analysis of the case where the colors are symmetry related, e.g. spin up and down in magnetic systems. In spite of their efficiency in counting SICs, neither Pólya's nor De Bruijn's theory helps in solving the difficult problem of identifying the complete list of SICs. Representative SICs are obtained by adopting an orderly generation approach, based on lexicographic ordering, which offers the advantage of avoiding the (computationally expensive) analysis and storage of all the possible configurations. When the number of colors increases, this strategy can be combined with the surjective resolution principle, which permits the efficient generation of SICs of a problem in |R| colors starting from the ones obtained for the (|R| - 1)-colors case. The whole scheme is documented by means of three examples: the abstract case of a square with C(4v) symmetry and the real cases of the garnet and olivine mineral families.

Ab initio modeling of trititanate nanotubes.

Nanotubes with H2Ti3O7 stoichiometry have been studied by ab initio methods in the diameter range of ∼20-60 Å. The distribution of the protons on the unit cell surface and its effect on the structure and electronic properties of the titanium dioxide backbone has been investigated. The preferred protonation pattern was found to be the function of the diameter of the nanotube and the convexity of the surface. Trititanate tubes have been found to be more stable than lepidocrocite nanotubes at all diameters. At diameters<50 Å the dissociation of water from the internal surface and the emergence of a new tube wall structure was observed. Band gaps were found to change only slightly with the curvature, the maximum difference from the flat slab values was ±0.5 eV.

On the use of symmetry in the ab initio quantum mechanical simulation of nanotubes and related materials.

Nanotubes can be characterized by a very high point symmetry, comparable or even larger than the one of the most symmetric crystalline systems (cubic, 48 point symmetry operators). For example, N = 2n rototranslation symmetry operators connect the atoms of the (n,0) nanotubes. This symmetry is fully exploited in the CRYSTAL code. As a result, ab initio quantum mechanical large basis set calculations of carbon nanotubes containing more than 150 atoms in the unit cell become very cheap, because the irreducible part of the unit cell reduces to two atoms only. The nanotube symmetry is exploited at three levels in the present implementation. First, for the automatic generation of the nanotube structure (and then of the input file for the SCF calculation) starting from a two-dimensional structure (in the specific case, graphene). Second, the nanotube symmetry is used for the calculation of the mono- and bi-electronic integrals that enter into the Fock (Kohn-Sham) matrix definition. Only the irreducible wedge of the Fock matrix is computed, with a saving factor close to N. Finally, the symmetry is exploited for the diagonalization, where each irreducible representation is separately treated. When M atomic orbitals per carbon atom are used, the diagonalization computing time is close to Nt, where t is the time required for the diagonalization of each 2M x 2M matrix. The efficiency and accuracy of the computational scheme is documented.

Single-layered chrysotile nanotubes: A quantum mechanical ab initio simulation.

Chrysotile single-layered nanotubes, obtained by wrapping the Mg(3)Si(2)O(5)(OH)(4) lizardite monolayer along the (n,-n) hexagonal lattice vector, are simulated at the ab initio level by using an all electron 6-31G( *) basis set and the B3LYP functional for n varying from 14 to 24 (the nanotube radius R referred to the oxygen connecting the Mg and Si layers increases from 20 to 35 A). Because of the full exploitation of the helical symmetry, recently implemented in the CRYSTAL code, the computational cost for the full self-consistent field (SCF) and gradient calculation increases only by a factor of 2 and 1.2, respectively, when passing from the lizardite monolayer [18 atoms and 236 AOs (atomic orbitals) in the unit cell] to the (24, -24) tube (864 atoms and 11,328 AOs). The total energy of the tubes is always larger than that of the lizardite monolayer; the difference DeltaE decreases very rapidly with n; for the largest tube here considered (n=24) DeltaE is as small as 2.7 kJ/mol per formula unit (f.u.); extrapolating to larger n values, at about R=50 A, DeltaE becomes smaller than 1 kJ mol f.u. Very large energy gains are observed for small n values during optimization after rolling, mainly due to the rotation of the SiO(4) tetrahedra that are in the inner part of the cylinder ("normal rolling"); such a rigid rotation accounts for about 85% of the overall relaxation energy. "Inverse rolling" tubes (SiO(4) on the external wall of the tube) are shown to be less stable than the corresponding "normal" tubes.